The bignum Library

This library is an adaptation of the “miracl” bignum library published by Mike Scott in 1989. Since then, Mike has developed his library further. The latest version can be downloaded from his web site

ftp.compapp.dcu.ie/pub/crypto/miracl.zip

I had worked a lot in the 1989 version, rewriting the essential parts of it in pure assembly. It runs quite fast now, and I am convinced it is a very useful package.

The interface this library presents is described in “bignums.h”. Basically, most operators (addition, subtraction, etc.) are overloaded there, so that you can use the library using a practical infix notation (you write c = a + b instead of quadsum(a,b,c) to make an addition).

Practical considerations:

Note that the bignum library is designed to work with the garbage collector. Most operators that need to, for instance “+”, “-“ etc., return a new bignum, that will be garbage collected automatically if it is no longer needed. This can be costly, and if you prefer trading speed for notational convenience, you can use the normal C interface and instead of writing:

c = (a + b)/(a – b);//creates temporary objects for a+b and a-b

you write:

c = newBignum(0); //allocate result

quadsum(a,b,tmp); // tmp = a + b

quadsub(a,b,c); // c = a - b

quaddiv(tmp,c,c); // c = tmp/c = (a+b)/(a-b)

In this latter case, you create only a single temporary that can later be reused. In the former case, you create two. This is absolutely unnoticeable in MOST cases. The convenience of the infix notation largely compensates the occasional stops to collect unneeded structures. Of course, you can trade execution speed for development speed.

Note that the public interface of the bignum library is an opaque data structure:

typedef struct _mp {

void *mp;

} pBignum;

This allows for changes in the bignum representation to be made transparently to your code. You could plug in a completely different bignum library, or get the latest version of Mike Scott and plug it in, without any changes to the bignum client code. Note that this is just a glorified pointer. Before using it, you must initialize it to a new bignum, either as a result of an operation or as the result of newBignum().

The library is presented as a DLL. When loaded, the DLL will automatically initialize the bignum system to a precision of 300 bits. You can change this using the BignumPrecision() API.

Example:

This example prints the powers of two from 127 to zero.

#include <bignums.h>

int main(void)

{

pBignum a,b,c;

char buffer[4096];

c = 2; // Note the overloaded assignment

for (int i = 127;i> 0;i--) {

b = newBignum(i);

a = quadpow(c,b);

quadformat(a,buffer);

printf("[%3d] %s\n",i,buffer);

}

return 0;

}

Output:

[127] 170141183460469231731687303715884105728

[126] 85070591730234615865843651857942052864

[125] 42535295865117307932921825928971026432

[124] 21267647932558653966460912964485513216

[123] 10633823966279326983230456482242756608

[122] 5316911983139663491615228241121378304

[121] 2658455991569831745807614120560689152

[120] 1329227995784915872903807060280344576

[119] 664613997892457936451903530140172288

…

Interface overview

Function	Description
Initialization
int BignumPrecision(int);	Sets the precision to a given value. Each big number is represented as a vector of integers. For each integer, the precision increases of 30 bits. The argument is the number of integers to reserve for each big number. For instance, if you request a precision of 30, 30 integers of 30 bits precision will be used, i.e. 900 bits per big number.
Creation
pBignum newBignum(void)	Returns a new big number initialized to zero
pBignum newBignum(int);	Returns a new bignum initialized to the given integer.
pBignum newBignum(double);	Returns a new big number initialized to the given double precision number.
pBignum newBignum(long long);	Returns a new big number initialized to the given long long integer.
pBignum newBignum(char *);	Returns a new big number initialized to the converted value of the given string.
Conversions
pBignum atoquad(char *);	Converts the given string to a newly allocated big number.
pBignum longToquad(long); pBignum doubleToquad(double); pBignum longlongToquad(long long);	Converts the given 32 bit integer,double precision, or long long number to a newly allocated big number.
void long2quad(long,pBignum); void double2quad(double,pBignum);	Converts the given 32 bit integer or double precision number into a big number. Sotrage is not allocated, the result is stored in the given bignum, that must have been previously allocated.
double quad2double(pBignum); long long quadToll(pBignum);	Returns a double precision or a long long number from a big number.
Formatting
void quadexpformat(pBignum x, unsigned char *outbuf, int width, int decimals);	Formats a given big number “x” as text in the buffer pointed to by “outbuf”, with the given width and decimals.
void quadtoa(pBignum,unsigned char *outbuf);	Formats a given big number as text in the given output buffer.
Comparisons
int operator<(pBignum a,pBignum b); int operator<(pBignum,long); int operator<(pBignum,double); int operator<(long,pBignum); int operator<(double,pBignum);	(a<b)?1:0
int operator<=(pBignum,pBignum); int operator<=(pBignum,long); int operator<=(pBignum,double); int operator<=(long,pBignum); int operator<=(double,pBignum);	(a<=b)?1:0
int operator==(pBignum,pBignum); int operator==(pBignum,long); int operator==(pBignum,double); int operator==(long,pBignum); int operator==(double,pBignum);	(a==b)?1:0;
int operator>=(pBignum,pBignum); int operator>=(pBignum,long); int operator>=(pBignum,double); int operator>=(long,pBignum); int operator>=(double,pBignum);	(a>=b)?1:0;
int operator>(pBignum,pBignum); int operator>(pBignum,long); int operator>(pBignum,double); int operator>(long,pBignum); int operator>(double,pBignum);	(a>b)?1:0;
int operator!=(pBignum,pBignum); int operator!=(pBignum,long); int operator!=(pBignum,double); int operator!=(long,pBignum); int operator!=(double,pBignum);	(a!=b)?1:0;
pBignum operator=(pBignum &,int); pBignum operator=(double &,pBignum); pBignum operator=(pBignum &,double);	Assignments from int, double or to double.
int quadcmp(pBignum a,pBignum b)	Returns –1 if a < b, 0 if a == b, and 1 if a > b
Addition
operator+(pBignum,pBignum); pBignum operator+(pBignum,long); … other conversions	Returns newly allocated big number with the result of the addition.
operator+=(pBignum,pBignum);	Returns the addition of the first operand and the second. No new storage, the result is stored in the first operator.
Sustraction
operator-(pBignum,pBignum);	Returns newly allocated big number with the result of the sustraction.
operator-=(pBignum,pBignum);	Returns a pointer to the first argument. The substraction is stored in the first argument, no new storage is used.
Multiplication
operator *(pBignum,pBignum);	Returns newly allocated big number with the result of the multiplication.
operator *=(pBignum,pBignum);	Stores the multiplication result in its first argument. Returns first argument.
Division
operator /(pBignum,pBignum);	Returns newly allocated big number with the result of the division.
operator /=(pBignum,pBignum);	Stores the division result in its first argument. Returns first argument..
Math functions
pBignum quadisqrt(pBignum);	Returns a newly allocated big number with the integer square root of its argument.
pBignum quadsqrt(pBignum);	Returns a newly allocated big number with the square root of its argument.
pBignum logb2(pBignum);	Returns a newly allocated big number with the log base 2 of its argument.
pBignum quadln(pBignum);	Returns the natural logarithm of its arguments in a newly allocated big number.
pBignum quadpow(pBignum,pBignum) pBignum quadpow(pBignum,int);	Returns a newly allocated big number with the result of elevating the first argument to the second argument power, i.e. result = first ^Second
pBignum quadexp(pBignum)	Returns newly allocated big number containing e raised to the first argument.
pBignum quadrand(pBignum)	Returns a random number between 0 and the given big number
pBignum quadpgcd(pBignum x, pBignum y);	Returns the gcd of the two big numbers
pBignum quadround(pBignum x,int n);	Rounds the given number to n places. Returns a new number, the input is not modified.
pBignum quadratio(int numerator, int denominator);	Returns a new big number built as a ratio from the given integers.
int quadintegerp(pBignum x);	Returns 1 if the given big number is represented as an exact integer, or 0 if it is represented as a ratio.
pBignum quadfloor(pBignum x);	Returns the largest integer less than x
pBignum quadtrunc(pBignum x);	Return the nearest integer to x but not larger.
pBignum fact(int n);	Returns n!
Trig Note: All arguments must be in radians.
pBignum quadsin(pBignum);	Sinus
pBignum quadcos(pBignum);	Cosinus
pBignum quadacos(pBignum);	Arc cosinus
pBignum quadasin(pBignum);	Arc sinus
pBignum quadtan(pBignum);	Tangent
pBignum quadatan(pBignum);	Arc tangent
pBignum quadtanh(pBignum);	Hyperbolic tangent
pBignum quadatanh(pBignum);	Hyperbolic arc tangent
pBignum quadsinh(pBignum);	Hyperbolic sinus
pBignum quadasinh(pBignum);	Hyperbolic arc sinus
pBignum quadcosh(pBignum);	Hyperbolic cosinus
pBignum quadacosh(pBignum);	Hyperbolic arc cosinus

Debugger support.

Bignums are recognized in the latest versions of wedit and the debugger shows them as numbers.

Important considerations when writing bignum code

· A pBignum is actually a pointer to a number. This is more efficient, but can lead to mistakes. When you write:

pBignum a;

// Some code here

pBignum b = a;

After this assignment both a and b point to the same number! If you modify b, a will get modified too. To avoid this problem you should write:

pBignum b = quadcopy(a);

This forces b to point to a different area of memory, i.e. a copy of a. This semantics are obvious for double precision numbers or for integers, but here they must be explicitly enforced by the programmer.

· Avoid floating point numbers when working with bignums if possible. For instance, do not write:

pBignum a = 0.5; // WRONG

You should write:

pBignum a = quadratio(1,2); //OK

The later expression stores exactly 1/2. The former stores an approximation to 1/2 that will be good in the first 16 digits only.

Bignum implementation

Bignums are stored in two ways:

· As a large integer: The first word (32 bits) contains the length and the sign, followed by the data, 30 bits for each 32 bit word.

· As a ratio of two integers, represented as a numerator and denominator. This allows for a floating point (or maybe floating slash) representation.

In both cases, the library returns a pointer to these opaque data. No direct access is provided, so the underlying library or DLL can be changed as needed, without affecting user code, as long as the interface of the dll is the same.

Here is an excerpt from the "Miracl" reference manual that explains the basic implementation of the library.

Floating-Slash numbers

The straightforward way to represent rational numbers is as reduced fractions, as a numerator and denominator with all common factors cancelled out. These numbers can then be added, subtracted, multiplied and divided in the obvious way and the result reduced by dividing both numerator and denominator by their Greatest Common Divisor. An efficient GCD subroutine, using Lehmers modification of the classical euclidean algorithm for multiprecision numbers [Knuth81], is included in the MIRACL package.

An alternative way to represent rationals would be as a finite continued fraction [Knuth81]. Every rational number p/q can be written as

or more elegantly as p/q = [a₀/a₁/a₂/..../a_n] where the a_i are positive integers, usually quite small.

For example

Note that the a_i elements of the above continued fraction representation are easily found as the quotients generated as a by-product when the euclidean GCD algorithm is applied to p and q.

As we are committed to fixed length representation of rationals, a problem arises when the result of some operation exceeds this fixed length. There is a necessity for some scheme of truncation, or rounding. While there is no obvious way to truncate a large fraction, it is a simple matter to truncate the continued fraction representation. The resulting, smaller, fraction is called a best rational approximation, or a convergent, to the original fraction.

Consider truncating 277/642 = [0/2/3/6/1/3/3]. Simply drop the last element from the CF representation, giving [0/2/3/6/1/3] = 85/197, which is a very close approximation to 277/642 (error = 0.0018%). Chopping more terms from the CF expansion gives the successive convergents as 22/51, 19/44, 3/7, 1/2, 0/1. As the fractions get smaller, the error increases. Obviously the truncation rule for a computer implementation should be to choose the biggest convergent that fits the computer representation.

The type of rounding described above is also called ‘Mediant rounding’. If p/q and r/s are two neighbouring representable slash numbers astride a gap, then their mediant is the unrepresentable (p+r)/(q+s). All larger fractions between p/q and the mediant will round to p/q, and those between r/s and the mediant will round to r/s. The mediant itself rounds to the ‘simpler’ of p/q and r/s.

This is theoretically a very good way to round, much better than the rather arbitrary and base-dependent methods used in floating-point arithmetic, and is the method used here. The full theoretical basis of floating-slash arithmetic is described in detail by Matula & Kornerup [Matula85]. It should be noted that our flash representation is in fact a cross between the fixed- and floating-slash systems analysed by Matula & Kornerup, as our slash can only float between words, and not between bits. However the characteristics of the flash data-type will tend to those of floating-slash, as the precision is increased.

The MIRACL routine round implements mediant rounding. If the result of an arithmetic operation is the fraction p/q, then the euclidean GCD algorithm is applied as before to p and q. However this time the objective is not to use the algorithm to calculate the GCD per se, but to use its quotients to build successive convergents to p/q. This process is stopped when the next convergent is too large to fit the flash representation. The complete algorithm is given below (Kornerup & Matula [Korn83])

Given p³0 and q³1

b_-2=p x_-2=0 y_-2=1

b_-1=q x_-1=1 y_-1=0

Now for i=0,1,..... and for b_i-1>0, find the quotient a_i and remainder b_i when b_i-2 is divided by b_i-1, such that

b_i = -a_i.b_i-1 + b_i-2

Then calculate

x_i = a_i.x_i-1 + x_i-2

y_i = a_i.y_i-1 + y_i-2

Stop when is to big to fit the flash representation, and take as the rounded result.

If applied to 277/642, this process will give the same sequence of convergents as stated earlier.

Since this rounding procedure must be applied to the result of each arithmetic operation, and since it is potentially rather slow, a lot of effort has been made to optimise its implementation. Lehmer's idea of operating only with the most significant piece of each number for as long as possible [Knuth81] is used, so that for most of the iterations only single-precision arithmetic is needed. Special care is taken to avoid the rounded result overshooting the limits of the flash representation [Scott89a]. The application of the basic arithmetic routines to the calculation of elementary functions such as log(x), exp(x), sin(x), cos(x), tan(x) etc., uses the fast algorithms described by Brent [Brent76].

A disadvantage of using a flash type of variable to approximate real arithmetic is the non-uniformity in gap-size between representable values (Matula & Kornerup [Matula85]).

To illustrate this consider a floating-slash system which is constrained to have the product of numerator and denominator less than 256. Observe that the first representable fraction less than 1/1 in such a system is 15/16, a gap of 1/16. The next fraction larger than 0/1 is 1/255, a gap of 1/255. In general, for a k-bit floating-slash system, the gap size varies from smaller than 2^-k to a worst case 2^-k/2. In practise this means that a real value that falls into one of the larger gaps, will be represented by a fraction which will be accurate to only half its usual precision. Fortunately such large gaps are rare, and increasingly so for higher precision, occurring only near simple fractions. However it does mean that real results can only be completely trusted to half the given decimal places. A partial solution to this problem would be to represent rationals directly as continued fractions. This gives a much better uniformity of gap-size (Kornerup & Matula [Korn85]), but would be very difficult to implement using a high level language.

Arithmetic on flash data-types is undoubtedly slower than on an equivalent sized multiprecision floating-point type (e.g. [Brent78]). The advantages of the flash approach are its ability to exactly represent rational numbers, and do exact arithmetic on them. Even when rounding is needed, the result often works out correctly, due to the tendency of mediant-rounding to prefer a simple fraction over a complex one. For example the roots program (Chapter 8) when asked to find the square root of 2 and then square the result, comes back with the exact answer of 2, despite much internal rounding.